テーブルデータを行列とみなし、各種の線形代数手法を適用することで、データの特徴の可視化、潜在因子の抽出、容量圧縮などができる。その結果は分析者個人に依らない。ならば1変数、2変数のデータから基本的な統計量、分布、散布図を示すことがマナーであるように、多変数のテーブルデータから、その相関行列、主成分分析または特異値分解を示すことも基本的な記述統計ではないかと考える。データ分析の前段階としてそれらの十分な記述統計を準備することで、データから仮説検証するときも、予測モデルをつくるときも、(経験、勘でなく)その分析方法をより科学的な根拠から説明することができるだろう。
ここではテーブルデータの特異値分解の基本についてまとめる。特異値分解は画像分析、音声分析、リコメンデーションなど実社会のビッグデータに対して最も適用されている分析手法である。今回は小さなデータを用いた説明になるが、小さな行列で成立する方法をそのまま何万次元のデータにも適用できることが線形代数手法の長所でもある。なお特異値分解を含め、行列分解の概要については別の小文でもまとめている。
(https://kng-mtd.github.io/test0/matrixDecomposition)
grViz("digraph {
graph [label='観測値と因子', labelloc=t,
layout = dot, rankdir = TB,
nodesep=0.5,ranksep=1]
node [shape=rectangle]
f1 [label = '体力',shape=ellipse]
f2 [label = '知力',shape=ellipse]
f3 [label = '運',shape=ellipse]
v1 [label = 'マージャン']
v2 [label = '年俸']
v3 [label = '100m走']
v4 [label = '遠投']
f1->v1
f1->v2
f1->v3
f1->v4
f2->v1
f2->v2
f2->v3
f2->v4
f3->v1
f3->v2
f3->v3
f3->v4
}")
| 選手 | 100m走 | 遠投 | マージャン | 年棒 |
しおみ ながおか やまだ むらかみ オスナ サンタナ なかむら やまざき いしかわ |
6.0 6.1 6.5 7.0 7.7 7.5 7.5 6.5 7.7 |
100 80 70 70 60 50 80 100 70 |
3 4 6 4 8 6 8 4 3 |
10,000 5,000 20,000 20,000 15,000 15,000 12,000 3,000 7,000 |
たとえば、スポーツ選手を比較しようとするとき、まずいろいろな観測値を集める。それである程度の各選手の特徴を把握できる。しかし観測値は種類が多く(いくらでも集められる)、また観測値は多様な因子の組み合わせで決まるので、ある観測値が大きいと別の観測値も大きい、小さいなどの複雑な関係がある。分析者により、どの観測値を選び、どのように組み合わせるかに差、クセがでてくる、「いや、オレはこの数字とこの数字を見るよ」。
もっと一般的に、体系的に対象の特徴抽出、比較、分類するには、複雑な関係にある観測値ではなく、観測値の元となる、かつ互いに関係を持たない因子を見つけるといい。因子間に関係がなければ、単純に個々の因子を大きい順に並べて、わかりやすく比較することができる。各因子が互いに関係を持たないとき(無相関のとき)、それらを「互いに直交する」という。主成分分析、特異値分解はそのような直交する因子をみつけて、対象を比較、分類、可視化する手法である(データ特有のハイパーパラメータなども不要)。
特異値分解は、データ行列から互いに直交する因子を抽出し、各因子の重要性、各対象の特徴、各変数の特徴を同時に可視化するものである。対象を行、その変数を列とするデータ行列を、対象の特徴を表す直交する因子(左特異ベクトル)、変数の特徴を表す直交する因子(右特異ベクトル)、因子の重要性を表す特異値行列に分解するもの。
主成分分析は、データがもつ変数の相関行列から互いに直交する因子(主成分)を抽出し、各因子の重要性、各対象の因子の大きさ(主成分スコア)を可視化する。
変数間の相関行列を固有値分解し、固有ベクトル=主成分と、固有値=主成分の重要性、各対象がもつ主成分方向の大きさ(主成分スコア)を計算する。数学的には、互いに直行する分散の大きい方向をみつけることは、相関行列の固有値分解問題であり、対称行列である相関行列がもつ固有ベクトルは直交することが示されている。
前処理により特異値分解、主成分分析は対象、変数の結果は同じになる(どちらでもいい)。計算が早い特異値分解、結果を理解しやすい主成分分析あたりで使い分けられている。
例として3つの島に生息する、3種のペンギンの性別、体重、羽の長さ、嘴の長さ、嘴の幅の3年間調査であるpalmerpenguinsデータを使う。
tb0=penguins
glimpse(tb0)
## Rows: 344
## Columns: 8
## $ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex <fct> male, female, female, NA, female, male, female, male…
## $ year <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…
summary(tb0)
## species island bill_length_mm bill_depth_mm
## Adelie :152 Biscoe :168 Min. :32.1 Min. :13.1
## Chinstrap: 68 Dream :124 1st Qu.:39.2 1st Qu.:15.6
## Gentoo :124 Torgersen: 52 Median :44.5 Median :17.3
## Mean :43.9 Mean :17.1
## 3rd Qu.:48.5 3rd Qu.:18.7
## Max. :59.6 Max. :21.5
## NA's :2 NA's :2
## flipper_length_mm body_mass_g sex year
## Min. :172 Min. :2700 female:165 Min. :2007
## 1st Qu.:190 1st Qu.:3550 male :168 1st Qu.:2007
## Median :197 Median :4050 NA's : 11 Median :2008
## Mean :201 Mean :4202 Mean :2008
## 3rd Qu.:213 3rd Qu.:4750 3rd Qu.:2009
## Max. :231 Max. :6300 Max. :2009
## NA's :2 NA's :2
#datatable(tb0,options=list(scrollX=500))
欠損値のある対象を除く。年度を量的変数から質的変数に型変換する。
tb0 %<>%
drop_na() %>%
mutate(year=as.factor(year))
量的変数の正規化
特異値分解でも主成分分析と同様、変数の[0,1]正規化や標準正規化をする、変数の単位が異なるときは必須。一部の変数を正規化することが障害になるかどうか、特異値分解前の変数の記述統計により検討し、その変数を除くことで対処する。もし正規化しない場合、特異値分解の結果の解釈を難しくする。
tbn=tb0 %>%
select(is.numeric) # just numeric variable
nn=names(tbn)
l=c()
u=c()
for(i in nn){
l[i]=min(tbn[i])
u[i]=max(tbn[i])
tbn[i]=scale(tbn[i],center=l[i],scale=u[i]-l[i])
}
質的変数のダミー変数化
質的変数を含むデータに対しても線形代数手法を適用するために、2値の質的変数は{0,1}のダミー変数、3水準以上の質的変数は複数の{0,1}のダミー変数に変換する。また質的変数が順序をもつとき、整数値の連続変数に変換することでほぼ対処できる。
factor2ind=function(x, baseline){
xname=deparse(substitute(x))
n=length(x)
x=as.factor(x)
if(!missing(baseline)) x=relevel(x, baseline)
X=matrix(0L, n, length(levels(x)))
X[(1:n) + n*(unclass(x)-1)]=1L
X[is.na(x),]=NA
dimnames(X)=list(names(x), paste(xname, levels(x), sep = ":"))
return(X[,-1,drop=FALSE])
}
tbc0=tb0 %>%
select(!is.numeric) # just categorical variable
nc=names(tbc0)
tbc=tibble(rep(NA,nrow(tbc0)))
for(i in nc){
btb=as_tibble(factor2ind(tbc0[[i]]))
names(btb)=names(btb) %>%
str_replace('tbc0\\[\\[i\\]\\]',i)
tbc=tbc %>%
bind_cols(btb)
}
tbc=tbc[,-1]
すべての変数から成る相関行列をみて、他の変数との相関が微小な連続変数を除く(当然、質的変数の別水準を表すダミー変数と逆相関があるだけのダミー変数も除く)。そのような変数は特異値分解の結果に影響しない。
tb=bind_cols(tbn,tbc)
par(mfrow=c(1,1))
corrplot(cor(tb))
tb=tb %>%
select(-'year:2008',-'year:2009')
tbc=tbc %>%
select(-'year:2008',-'year:2009')
特異値分解する、前処理後のデータ行列
glimpse(tb)
## Rows: 333
## Columns: 9
## $ bill_length_mm <dbl> 0.2545, 0.2691, 0.2982, 0.1673, 0.2618, 0.2473, 0.…
## $ bill_depth_mm <dbl> 0.667, 0.512, 0.583, 0.738, 0.893, 0.560, 0.774, 0…
## $ flipper_length_mm <dbl> 0.1525, 0.2373, 0.3898, 0.3559, 0.3051, 0.1525, 0.…
## $ body_mass_g <dbl> 0.292, 0.306, 0.153, 0.208, 0.264, 0.257, 0.549, 0…
## $ `species:Chinstrap` <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ `species:Gentoo` <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ `island:Dream` <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ `island:Torgersen` <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0,…
## $ `sex:male` <int> 1, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1,…
summary(tb)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.269 1st Qu.:0.298 1st Qu.:0.305 1st Qu.:0.236
## Median :0.451 Median :0.500 Median :0.424 Median :0.375
## Mean :0.432 Mean :0.484 Mean :0.491 Mean :0.419
## 3rd Qu.:0.600 3rd Qu.:0.667 3rd Qu.:0.695 3rd Qu.:0.576
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
## species:Chinstrap species:Gentoo island:Dream island:Torgersen
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :0.000 Median :0.000 Median :0.000
## Mean :0.204 Mean :0.357 Mean :0.369 Mean :0.141
## 3rd Qu.:0.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:0.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
## sex:male
## Min. :0.000
## 1st Qu.:0.000
## Median :1.000
## Mean :0.505
## 3rd Qu.:1.000
## Max. :1.000
特異値分解では、データ行列(行が各対象、列は変数)の元の変数と同じ数(本データは9コ)の特異値が計算される。 適当な「次数」=因子数を決めて、特異値分解を打ち切る。対象数×次数の左特異ベクトル行列、次数次元の正方行列である特異値行列、変数の数×次数の右特異ベクトル行列ができる。
たとえば次数を4としたとき、左特異ベクトル行列は対象を4コの数値(=因子の大きさ)で表したもの。特異値行列は各因子が対象の特徴をどれだけ表すかの重み。右特異ベクトル行列は各因子が元の各変数とどれだけ相関するかを示す重みであり、すなわち元の変数の特徴を4コの数値で表したものとなる。
k=4
ku=k
kv=k
s=svd(tb,ku,kv) #u:n x ku, d:ku x kv, v:m x kv
#s
#s$d #singular value
u=s$u #left singular matrix
head(u,10)
## [,1] [,2] [,3] [,4]
## [1,] -0.0458 0.00767 0.1258 -0.0539
## [2,] -0.0243 0.00181 0.0571 -0.1394
## [3,] -0.0262 0.00308 0.0569 -0.1454
## [4,] -0.0269 0.00615 0.0638 -0.1493
## [5,] -0.0517 0.01009 0.1297 -0.0669
## [6,] -0.0225 0.00498 0.0603 -0.1386
## [7,] -0.0555 0.00179 0.1239 -0.0661
## [8,] -0.0219 0.00592 0.0585 -0.1383
## [9,] -0.0533 0.01063 0.1316 -0.0701
## [10,] -0.0554 0.00613 0.1314 -0.0717
d=diag(s$d[1:k])
d
## [,1] [,2] [,3] [,4]
## [1,] 23.1 0.0 0.00 0.00
## [2,] 0.0 13.1 0.00 0.00
## [3,] 0.0 0.0 8.94 0.00
## [4,] 0.0 0.0 0.00 6.74
v=s$v #right singular matrix
head(v,10)
## [,1] [,2] [,3] [,4]
## [1,] -0.3620 -0.0313 -0.1654 -0.0997
## [2,] -0.3633 0.2616 0.2438 -0.2756
## [3,] -0.4026 -0.2163 -0.1358 -0.1665
## [4,] -0.3501 -0.2091 -0.0413 -0.0497
## [5,] -0.1816 0.3939 -0.3387 -0.0689
## [6,] -0.3078 -0.5738 -0.3605 -0.0478
## [7,] -0.2911 0.5961 -0.3272 0.0387
## [8,] -0.0762 0.0134 0.4749 -0.7163
## [9,] -0.4838 0.0147 0.5620 0.6019
元の行列は左特異ベクトル行列、特異値行列、右特異ベクトル行列の積で近似される(近似行列)。特異値分解による近似行列は、あらゆる行列分解のうち、元の行列に最も近似する(行列要素間の差のフロベニウスノルム最小)(Eckart-Youngの定理)。左右の特異ベクトル行列、特異値行列の要素数の合計が、元のデータ行列の要素数よりも少なくなれば圧縮できるが、どれだけ元の行列と近似行列が近いのかを確認するには、以下のような方法があるだろう。
行列全体の数値をヒートマップで(直感的に)比べる
量的変数の各要素の差をみる(平均二乗誤差 MSE、二乗誤差の分散 V[SE])
近似行列の各量的変数の基本統計量をみる
各質的変数のクロス表をみる(近似行列の数値は0.5を境界に0,1に分ける)
各量的変数の分布をボックスプロット、ヒストグラムで(直感的に)比べる
各量的変数の平均の差をt検定する
相関行列を比べる
散布図行列を(直感的に)比べる
row.names(v)=names(tb)
tb_svd0=as_tibble(u %*% d %*% t(v))
tbn_svd=tb_svd0 %>%
select(names(tbn))
tbc_svd=tb_svd0 %>%
select(names(tbc))
tbc_svd[tbc_svd<.5]=0
tbc_svd[tbc_svd>.5]=1
tb_svd=bind_cols(tbn_svd,tbc_svd)
heatmap(as.matrix(tb),Rowv=NA,Colv=NA,scale='none',margins=c(10,3),main='origin')
#heatmap(u %*% d %*% t(v),Rowv=NA,Colv=NA,scale='none',margins=c(10,3))
heatmap(as.matrix(tb_svd),Rowv=NA,Colv=NA,scale='none',margins=c(10,3),main='SVD4')
cat('MSE: ',mean((tbn-tbn_svd)**2 |> unlist()))
## MSE: 0.00859
cat('V[SE]: ',var((tbn-tbn_svd)**2 |> unlist()))
## V[SE]: 0.000335
summary(tb_svd)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## Min. :0.108 Min. :0.039 Min. :0.107 Min. :0.051
## 1st Qu.:0.303 1st Qu.:0.301 1st Qu.:0.292 1st Qu.:0.217
## Median :0.478 Median :0.452 Median :0.417 Median :0.383
## Mean :0.439 Mean :0.462 Mean :0.489 Mean :0.415
## 3rd Qu.:0.576 3rd Qu.:0.627 3rd Qu.:0.705 3rd Qu.:0.563
## Max. :0.759 Max. :0.957 Max. :0.996 Max. :0.894
## species:Chinstrap species:Gentoo island:Dream island:Torgersen
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :0.000 Median :0.000 Median :0.000
## Mean :0.216 Mean :0.357 Mean :0.369 Mean :0.141
## 3rd Qu.:0.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:0.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
## sex:male
## Min. :0.000
## 1st Qu.:0.000
## Median :1.000
## Mean :0.505
## 3rd Qu.:1.000
## Max. :1.000
for(i in names(tbc)){
cat('\n',i)
table(tbc[[i]],tbc_svd[[i]]) %>% print()
#wilcox.test(tb[[i]],tb_svd[[i]]) %>% print()
}
##
## species:Chinstrap
## 0 1
## 0 261 4
## 1 0 68
##
## species:Gentoo
## 0 1
## 0 214 0
## 1 0 119
##
## island:Dream
## 0 1
## 0 210 0
## 1 0 123
##
## island:Torgersen
## 0 1
## 0 286 0
## 1 0 47
##
## sex:male
## 0 1
## 0 165 0
## 1 0 168
par(mfrow=c(1,1))
boxplot(tbn,main='origin')
boxplot(tbn_svd,main='SVD4')
par(mfrow=c(2,1))
for(i in names(tbn)){
hist(tb[[i]],main=i,xlim=c(0,1.5))
hist(tb_svd[[i]],main='',xlim=c(0,1.5))
#a=cor(tb[i],tb_svd[i])
#cat(i,' correlation:',a[[1]],'\n')
t.test(tbn[i],tbn_svd[i]) %>% print()
}
##
## Welch Two Sample t-test
##
## data: tbn[i] and tbn_svd[i]
## t = -0.4, df = 638, p-value = 0.7
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0339 0.0213
## sample estimates:
## mean of x mean of y
## 0.432 0.439
##
## Welch Two Sample t-test
##
## data: tbn[i] and tbn_svd[i]
## t = 1, df = 663, p-value = 0.2
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0140 0.0587
## sample estimates:
## mean of x mean of y
## 0.484 0.462
##
## Welch Two Sample t-test
##
## data: tbn[i] and tbn_svd[i]
## t = 0.1, df = 664, p-value = 0.9
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0334 0.0381
## sample estimates:
## mean of x mean of y
## 0.491 0.489
##
## Welch Two Sample t-test
##
## data: tbn[i] and tbn_svd[i]
## t = 0.2, df = 664, p-value = 0.8
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0301 0.0373
## sample estimates:
## mean of x mean of y
## 0.419 0.415
cat('origin')
## origin
cor(tbn)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm 1.000 -0.229 0.653 0.589
## bill_depth_mm -0.229 1.000 -0.578 -0.472
## flipper_length_mm 0.653 -0.578 1.000 0.873
## body_mass_g 0.589 -0.472 0.873 1.000
cat('SVD4')
## SVD4
cor(tbn_svd)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm 1.0000 -0.0389 0.817 0.763
## bill_depth_mm -0.0389 1.0000 -0.383 -0.304
## flipper_length_mm 0.8172 -0.3826 1.000 0.972
## body_mass_g 0.7630 -0.3038 0.972 1.000
plot(tbn,main='origin')
plot(tbn_svd,main='SVD4')
# psyck::r.test
# 1.H0 r12=0, n single group
# r.test(n,r12)
# 2,H0 r12=r34, n1,n2 not correspond groups
# r,test(n=n1,r12,n2=n2,r34)
# 3,H0 r12=r13, n correspond groups
# r.test(n,r12,r23,r13)
# 4.H0 r12=r34, n correspond groups
# r.test(n,r12,r34,r23,r13,r14,r24)
#rtest=r.test(nrow(tbn),cor(tbn),cor(tbn_svd))
#rtest$p
近似行列において、量的変数の平均、中央値は再現されている。ダミー変数の値{0,1}を、近似行列では{0,1}に近い数値でほぼ再現されている。分布の形はほぼ再現されるが、ノイズ、外れ値が除去されている。相関関係の大小がより強調されている。散布図から、複数のクラスの別々の相関関係をみることができる。
観測値を単純化した近似行列を、母集団の母数とすることができるかもしれない。観測値はそれにノイズが加わったものとするモデルが考えられる。そこで特異値分解の結果にノイズを加えた「観測値」をつくってみる。
近似行列の全要素に上下10%の一様乱数を加えたものを10コから成る、対象数を10倍にした拡張データをつくった。拡張データから元と同じ対象数をサンプリングをすることで「観測」をシミュレーションした。
expd0=tibble()
for(i in 1:10){ # 10times
tb1=as_tibble(u %*% d %*% t(v))*runif(nrow(tb)*ncol(tb),0.9,1.1) # ±10%
expd0=bind_rows(expd0,tb1)
}
expdn=expd0 %>%
select(names(tbn))
expdc=expd0 %>%
select(names(tbc))
expdc[expdc<.5]=0
expdc[expdc>.5]=1
expd=bind_cols(expdn,expdc)
#glimpse(expd)
datatable(expd,options=list(scrollX=500))
summary(expd)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## Min. :0.099 Min. :0.037 Min. :0.102 Min. :0.047
## 1st Qu.:0.296 1st Qu.:0.295 1st Qu.:0.291 1st Qu.:0.217
## Median :0.470 Median :0.452 Median :0.415 Median :0.381
## Mean :0.439 Mean :0.462 Mean :0.488 Mean :0.415
## 3rd Qu.:0.572 3rd Qu.:0.630 3rd Qu.:0.709 3rd Qu.:0.573
## Max. :0.820 Max. :1.043 Max. :1.060 Max. :0.966
## species:Chinstrap species:Gentoo island:Dream island:Torgersen
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.000
## 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000
## Median :0.000 Median :0.000 Median :0.000 Median :0.000
## Mean :0.231 Mean :0.357 Mean :0.369 Mean :0.141
## 3rd Qu.:0.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:0.000
## Max. :1.000 Max. :1.000 Max. :1.000 Max. :1.000
## sex:male
## Min. :0.000
## 1st Qu.:0.000
## Median :1.000
## Mean :0.505
## 3rd Qu.:1.000
## Max. :1.000
write_csv(expd,'penguins_expand.csv')
expd=read_csv('penguins_expand.csv')
tb1=sample_n(expd,333) #sample n=333
tb1n=tb1 %>% select(names(tbn))
par(mfrow=c(1,1))
boxplot(tbn,main='origin')
boxplot(tb1n,main='new obs.')
par(mfrow=c(2,1))
for(i in names(tbn)){
hist(tb[[i]],main=i,xlim=c(0,1.5))
hist(tb1[[i]],main='',xlim=c(0,1.5))
t.test(tb[i],tb1[i]) %>% print()
}
##
## Welch Two Sample t-test
##
## data: tb[i] and tb1[i]
## t = -0.5, df = 636, p-value = 0.6
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0345 0.0206
## sample estimates:
## mean of x mean of y
## 0.432 0.439
##
## Welch Two Sample t-test
##
## data: tb[i] and tb1[i]
## t = 0.6, df = 658, p-value = 0.6
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0264 0.0486
## sample estimates:
## mean of x mean of y
## 0.484 0.473
##
## Welch Two Sample t-test
##
## data: tb[i] and tb1[i]
## t = 0.8, df = 661, p-value = 0.4
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.02 0.05
## sample estimates:
## mean of x mean of y
## 0.491 0.476
##
## Welch Two Sample t-test
##
## data: tb[i] and tb1[i]
## t = 0.6, df = 660, p-value = 0.6
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.0233 0.0424
## sample estimates:
## mean of x mean of y
## 0.419 0.409
cat('origin')
## origin
cor(tbn)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm 1.000 -0.229 0.653 0.589
## bill_depth_mm -0.229 1.000 -0.578 -0.472
## flipper_length_mm 0.653 -0.578 1.000 0.873
## body_mass_g 0.589 -0.472 0.873 1.000
cat('new obs.')
## new obs.
cor(tb1n)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm 1.0000 0.0349 0.785 0.738
## bill_depth_mm 0.0349 1.0000 -0.353 -0.277
## flipper_length_mm 0.7848 -0.3526 1.000 0.951
## body_mass_g 0.7379 -0.2773 0.951 1.000
plot(tbn,main='origin')
plot(tb1n,main='new obs.')
library(imager)
img0=load.image('image1.jpg')
#plot(img0,axes=F)
img=grayscale(img0) #[0,1]
#plot(img,axes=F)
save.image(img,'image10.jpg')
mx=as.matrix(img)
img=as.cimg(mx)
plot(img,axes=F,main='origin')
w0=dim(mx)[1]
h0=dim(mx)[2]
cat('width: ',w0,' height: ',h0)
## width: 5472 height: 3648
img1=resize(img,w0/64,h0/64)
plot(img1,axes=F,main='1/64')
k=20
gc(reset=T)
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2554552 136 4124740 220 2554552 136
## Vcells 144330088 1101 267872290 2044 144330088 1101
#cat('\nSVD k:',k)
system.time({
s=svd(mx,k,k)
u=s$u
d=diag(s$d[1:k])
v=s$v
tb1=u %*% d %*% t(v)
})
## ユーザ システム 経過
## 92.51 0.64 195.02
img=as.cimg(tb1)
plot(img,axes=F,main=str_c('SVD ',k))
save.image(img,'image1svd20.jpg')
gc(reset=T)
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2555711 136 4124740 220 2555711 136
## Vcells 184443856 1407 329226953 2512 184443856 1407
#cat('\nrSVD k:',k)
system.time({
s=rsvd(mx,k)
u=s$u
d=diag(s$d[1:k])
v=s$v
tb1=u %*% d %*% t(v)
})
## ユーザ システム 経過
## 1.39 0.02 3.41
img=as.cimg(tb1)
plot(img,axes=F,main=str_c('rSVD ',k))
save.image(img,'image1rsvd20.jpg')
k=50
gc(reset=T)
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2557223 137 4124740 220 2557223 137
## Vcells 224366754 1712 395152343 3015 224366754 1712
#cat('\nSVD k:',k)
system.time({
s=svd(mx,k,k)
u=s$u
d=diag(s$d[1:k])
v=s$v
tb1=u %*% d %*% t(v)
})
## ユーザ システム 経過
## 72.92 0.38 184.48
img=as.cimg(tb1)
plot(img,axes=F,main=str_c('SVD ',k))
save.image(img,'image1svd50.jpg')
gc(reset=T)
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 2557437 137 4124740 220 2557437 137
## Vcells 244612716 1866 404516612 3086 244612716 1866
#cat('\nrSVD k:',k)
system.time({
s=rsvd(mx,k)
u=s$u
d=diag(s$d[1:k])
v=s$v
tb1=u %*% d %*% t(v)
})
## ユーザ システム 経過
## 2.52 0.10 6.03
img=as.cimg(tb1)
plot(img,axes=F,main=str_c('rSVD ',k))
save.image(img,'image1rsvd50.jpg')